import os
import pandas as pd
from datetime import datetime
import numpy as np
import plotly.graph_objects as go
import plotly.express as px
import plotly.io as pio
from gpx_file_reader import GPXFile

# this allows rendering of graphics when exported to HTML
pio.renderers.default = "plotly_mimetype+notebook"

Get the GPX file we need, and convert it to a gpxFile object using the parser.

GPX_FILE = '20210808_181500_amtrak_dc_to_nyc.gpx'
# print("abspath: ", abspath)
# dirname = os.path.dirname(abspath)
# print("dirname: ", dirname)
# os.chdir(dirname)
# files = os.listdir()
# print("files:\n", files)

gpx_file_path = os.path.abspath('20210808_181500_amtrak_dc_to_nyc.gpx')
gpxFile = GPXFile(gpx_file_path)

gpxFile.print_info()
gpxDF = gpxFile.get_gpx_dataframe()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 8684 entries, 0 to 8683
Data columns (total 4 columns):
 #   Column  Non-Null Count  Dtype         
---  ------  --------------  -----         
 0   lat     8684 non-null   float64       
 1   lon     8684 non-null   float64       
 2   ele     8684 non-null   float64       
 3   time    8684 non-null   datetime64[ns]
dtypes: datetime64[ns](1), float64(3)
memory usage: 271.5 KB
preview:
          lat        lon        ele                time
0  38.962250 -76.859985  24.902832 2021-08-08 22:15:01
1  38.962430 -76.859764  25.251038 2021-08-08 22:15:02
2  38.962654 -76.859510  24.751282 2021-08-08 22:15:03
3  38.962833 -76.859310  26.062988 2021-08-08 22:15:04
4  38.963060 -76.859055  29.826477 2021-08-08 22:15:05
5  38.963220 -76.858850  26.231506 2021-08-08 22:15:06
6  38.963436 -76.858610  22.677307 2021-08-08 22:15:07
7  38.963654 -76.858380  22.543457 2021-08-08 22:15:08
8  38.963852 -76.858154  22.253967 2021-08-08 22:15:09
9  38.964058 -76.857910  20.727356 2021-08-08 22:15:10

Plot elevation vs lat,lon

fig = go.Figure(data=go.Scattergeo(
        lon = gpxDF['lon'],
        lat = gpxDF['lat'],
        text = gpxDF['ele'],
        mode = 'markers',
        line = dict(width = 0.5,color = 'blue')
))

fig.update_layout(
        title = 'Amtrak train ride',
        geo_scope='usa',
        width=1000,
        height=800,
    )

fig.show()

Proof that much of the elevation was negative. (Inaccurate data)

fig = px.line_3d(gpxDF, x="lon", y="lat", z="ele")
fig.add_trace(go.Scatter3d(
    x = [gpxDF['lon'].iloc[0]], 
    y = [gpxDF['lat'].iloc[0]],
    z = [gpxDF['ele'].iloc[0]],
    name = 'DC'
))

fig.add_trace(go.Scatter3d(
    x = [gpxDF['lon'].iloc[-1]], 
    y = [gpxDF['lat'].iloc[-1]],
    z = [gpxDF['ele'].iloc[-1]],
    name = 'NYC'
))

fig.update_layout(
    title = "Elevation vs lat,lon",
    width=1000,
    height=800,
)
fig.show()

Plotting again, only taking into account positive elevation.

gpxDFElevationAbove0 = gpxDF[gpxDF['ele'] >= 0] 
fig = px.line_3d(gpxDFElevationAbove0, x="lon", y="lat", z="ele")
fig.add_trace(go.Scatter3d(
    x = [gpxDFElevationAbove0['lon'].iloc[0]], 
    y = [gpxDFElevationAbove0['lat'].iloc[0]],
    z = [gpxDFElevationAbove0['ele'].iloc[0]],
    name = 'DC'
))

fig.add_trace(go.Scatter3d(
    x = [gpxDFElevationAbove0['lon'].iloc[-1]], 
    y = [gpxDFElevationAbove0['lat'].iloc[-1]],
    z = [gpxDFElevationAbove0['ele'].iloc[-1]],
    name = 'NYC'
))

fig.update_layout(
    title = "Elevation vs lat,lon",
    width=1000,
    height=800,
)
fig.show()

# haversine formula: takes degrees
def dist_between(lat1, lon1, lat2, lon2):
    R = 6371e3 #metres
    phi1 = lat1 * np.pi/180
    phi2 = lat2 * np.pi/180
    deltaPhi = (lat2-lat1) * np.pi/180
    deltaLambda = (lon2-lon1) * np.pi/180
    
    a = np.sin(deltaPhi/2)**2 + np.cos(phi1) * np.cos(phi2) * (np.sin(deltaLambda/2)**2)
    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a))
    return R * c # meters

def test_dist_between():
    lat1 = gpxDF['lat'].iloc[0]
    lon1 = gpxDF['lon'].iloc[0]
    lat2 = gpxDF['lat'].iloc[-1]
    lon2 = gpxDF['lon'].iloc[-1]
    d = dist_between(lat1, lon1, lat2, lon2)
    print("dist (meters): ", d)
    dMiles = d/1609.34
    print("dist (miles): ", dMiles)
#test_dist_between()

# velocity calculation
gpxDF["deltaDistMeters"] = dist_between(gpxDF['lat'].shift(), gpxDF['lon'].shift(),
                                       gpxDF.loc[1:, 'lat'], gpxDF.loc[1:, 'lon'])
gpxDF["deltaTimeSeconds"] = (gpxDF.loc[1:, 'time'] - gpxDF['time'].shift()).apply(lambda row: row.total_seconds())
gpxDF["velocityMetersPerSecond"] = gpxDF["deltaDistMeters"] / gpxDF["deltaTimeSeconds"]
gpxDF["velocityMilesPerHour"] = gpxDF["velocityMetersPerSecond"] * (3600.0 / 1609.34)

gpxDFSpeedFiltered = gpxDF.loc[(gpxDF["velocityMilesPerHour"] >= 0) & (gpxDF["velocityMilesPerHour"] < 150)]
gpxDFSpeedFiltered.describe()

	lat	lon	ele	deltaDistMeters	deltaTimeSeconds	velocityMetersPerSecond	velocityMilesPerHour
count	8661.000000	8661.000000	8661.000000	8661.000000	8661.000000	8661.000000	8661.000000
mean	39.787519	-75.563328	-14.288275	35.543707	1.206905	34.525746	77.232086
std	0.432882	0.779921	10.459655	25.603769	5.413281	14.309065	32.008547
min	38.962430	-76.859764	-80.380710	0.222390	1.000000	0.131104	0.293272
25%	39.401394	-76.326996	-22.057251	22.881986	1.000000	22.592893	50.538988
50%	39.815760	-75.432526	-15.501739	37.802853	1.000000	37.569645	84.041111
75%	40.094870	-74.897545	-6.697693	46.452316	1.000000	46.197891	103.341996
max	40.584797	-74.303720	29.826477	1351.466293	303.000000	63.574951	142.213468

# Plot velocity

fig = px.line_3d(gpxDFSpeedFiltered, x="lon", y="lat", z="velocityMilesPerHour")
fig.add_trace(go.Scatter3d(
    x = [gpxDFSpeedFiltered['lon'].iloc[0]], 
    y = [gpxDFSpeedFiltered['lat'].iloc[0]],
    z = [gpxDFSpeedFiltered['velocityMilesPerHour'].iloc[0]],
    name = 'DC'
))

fig.add_trace(go.Scatter3d(
    x = [gpxDFSpeedFiltered['lon'].iloc[-1]], 
    y = [gpxDFSpeedFiltered['lat'].iloc[-1]],
    z = [gpxDFSpeedFiltered['velocityMilesPerHour'].iloc[-1]],
    name = 'NYC'
))

velocityColors = []
for i in range(0, len(gpxDFSpeedFiltered)):
    v = gpxDFSpeedFiltered['velocityMilesPerHour'].iloc[i]
    if v <= 50:
        color = 'red'
    elif 51 < v <= 100:
        color = 'yellow'
    elif 101 < v <= 150:
        color = 'green'
    else:
        color = 'black'
    velocityColors.append(color)

fig.add_trace(go.Scatter3d(
    x = gpxDFSpeedFiltered['lon'], 
    y = gpxDFSpeedFiltered['lat'],
    z = [0 for n in range(len(gpxDFSpeedFiltered))],
    name = 'Route',
    line=dict(color=velocityColors, width=0)
))
    
fig.update_layout(
    title = "Velocity vs lat,lon",
    width=1000,
    height=800,
)
fig.show()

# Distribution of velocities with histogram
fig = px.histogram(gpxDFSpeedFiltered, x="velocityMilesPerHour", nbins=20)
fig.show()

References:

GPX Schema documentation: http://www.topografix.com/GPX/1/1/#type_trksegType

GPX files with tag are corrupt:

https://github.com/NetTopologySuite/NetTopologySuite.IO.GPX/issues/34 This means that the GPX file I downloaded from the app actually does not conform fully to the schema.

Distance between 2 lat/lon coordinates: https://www.movable-type.co.uk/scripts/latlong.html

Vectorized haversine function: https://stackoverflow.com/questions/40452759/pandas-latitude-longitude-to-distance-between-successive-rows